#####################################
import warnings
warnings.filterwarnings("ignore")
import logging
logging.getLogger().setLevel(logging.CRITICAL)
PRINT=True
import h5py
import numpy as np
from numpy import linalg as LA
import matplotlib.pylab as plt
from matplotlib.colors import LogNorm
import illustris_python as il
import matplotlib.gridspec as gridspec
import matplotlib as mpl
from scipy.stats import median_abs_deviation as MAD # Median ABbsolute Deviation
%matplotlib inline
plt.style.use('dark_background')
plt.style.use('seaborn-bright')
COLOR = 'lightgray'
mpl.rcParams['text.color'] = COLOR # labels, titles, etc colors
mpl.rcParams['axes.labelcolor'] = COLOR # Numbers in axis colors
mpl.rcParams['xtick.color'] = COLOR # ticks colors
mpl.rcParams['ytick.color'] = COLOR # ticks colors
mpl.rc( 'axes',edgecolor = COLOR) # axis Colors
mpl.rc('font', size=15) # Font Size
mpl.rcParams['figure.figsize'] = [10, 10] # Size of figures
Data={}
LLL =[50,100,300]
VVV=[51.7,110.7,302.6]
h = 0.6774
POSITION= ['Fisrt','Second','Third']
# LLL =[50]
for i in LLL:
if PRINT:
print(f'Loading data of TNG{i}-1 ... ')
basePath = '/home/tnguser/sims.TNG/TNG'+str(i)+'-1/output/'
snapNum = 99 # z=0
fields = ['SubhaloFlag', 'SubhaloHalfmassRadType', 'SubhaloMassType', 'SubhaloPos', 'SubhaloVel',\
'SubhaloSFRinRad', 'SubhaloSFR','SubhaloStarMetallicity']
Subhalo = il.groupcat.loadSubhalos(basePath=basePath,snapNum=snapNum,fields=fields)
# SubhaloFlag = 0 means no cosmological origin
# stellar mass range log(M/h) = 9.0 to above
mask1= (Subhalo['SubhaloFlag']!=0)
idx = []
R200= []
M200= []
DM = Subhalo['SubhaloMassType'][:,1].copy()[mask1]
for PPP in range(3):
IDX=DM.argmax()
idx.append(IDX)
if PRINT:
print(f'In TNG{i}-1 the {POSITION[PPP]}\t most massive Halo ({IDX:^10}) has:\t {np.log10(DM[IDX]/h)+10:.3f} [Log(M_sun)] ')
DM[IDX]= 0
del DM
mask = np.where((Subhalo['SubhaloMassType'][:,4]>=0.003) & mask1)
globals()['mass_st_'+str(i)] = Subhalo['SubhaloMassType'][:,4][mask]/h # to save the mass in units of 10^10 solar masses
Data[str(i)]= dict( sh_ids = mask[0], # The subhalo "ID" is the same as the "index"
cm = Subhalo['SubhaloPos'][mask],
cv = Subhalo['SubhaloVel'][mask],
rh_type = Subhalo['SubhaloHalfmassRadType'][:,4][mask],
sfr_2rh = Subhalo['SubhaloSFRinRad'][mask], # in 2 * R_half
sfr_tot = Subhalo['SubhaloSFR'][mask], # in 2 * R_half
mh = Subhalo['SubhaloStarMetallicity'][mask],
MassiveH = np.array(idx))
# clear space
del Subhalo
if PRINT:
print('Done')
print()
Loading data of TNG50-1 ... In TNG50-1 the Fisrt most massive Halo ( 0 ) has: 14.225 [Log(M_sun)] In TNG50-1 the Second most massive Halo ( 63784 ) has: 13.897 [Log(M_sun)] In TNG50-1 the Third most massive Halo ( 96666 ) has: 13.757 [Log(M_sun)] Done Loading data of TNG100-1 ... In TNG100-1 the Fisrt most massive Halo ( 17066 ) has: 14.545 [Log(M_sun)] In TNG100-1 the Second most massive Halo ( 0 ) has: 14.537 [Log(M_sun)] In TNG100-1 the Third most massive Halo ( 31116 ) has: 14.464 [Log(M_sun)] Done Loading data of TNG300-1 ... In TNG300-1 the Fisrt most massive Halo ( 0 ) has: 15.205 [Log(M_sun)] In TNG300-1 the Second most massive Halo ( 11704 ) has: 15.047 [Log(M_sun)] In TNG300-1 the Third most massive Halo ( 17827 ) has: 14.974 [Log(M_sun)] Done
def mid(x):
return (x[1:]+x[:-1])/2
# From the graph
D1,D2 = np.genfromtxt('full_line.csv',delimiter=',',unpack=True)
l,u = np.log10(min(D1)),np.log10(max(D1)) # Range of masses desired
G1,G2 = np.genfromtxt('dots.dat',delimiter=',',unpack=True)
# From the NUMEX1
F1,F2,F3,F4,F5,F6 = np.genfromtxt('MyHMF.dat',delimiter='\t',unpack=True)
###################################################################
# plt.plot(D1,D2,label='Ilustris_TNG_paper',marker='.',color='k')
# plt.legend()
# plt.xscale('log')
# plt.yscale('log')
# plt.show(),plt.close()
# print(l,u)
fig= plt.figure(figsize=(18,6))
gs = gridspec.GridSpec( 1,5, # Number of axis y,x
height_ratios = [1], # relatives ratios of heigh
width_ratios = [1,0.1,1,0.1,1], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.15, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2],sharey=ax1,sharex=ax1)
ax3 = fig.add_subplot(gs[0,4],sharey=ax1,sharex=ax1)
AX=[ax1,ax2,ax3]
# AX[0].ticklabel_format(axis='y', scilimits=[-1, 1])
AX[0].set_ylim(1e-5,0.1)
# Ploting the masses
for j,i in enumerate(LLL):
VV = (VVV[j])**3
XX = np.log10(globals()['mass_st_'+str(i)])+10##np.append(log_mass_st_q,log_mass_st_sf)
XX = XX[(XX>l)&(XX<u)]
y,x = np.histogram(XX,bins=50)
w = x[1]-x[0]
AX[j].plot(10**mid(x),y/VV/w,label=f'TNG {i}-1',marker='s')
AX[j].plot(D1,D2,label='Ilustris_TNG_paper',marker='.',color='gray')
AX[j].scatter(G1,G2,label='GAMA-II (Wright et al. 2017)',marker='o',edgecolor='r',facecolor='none')
AX[j].legend(loc='upper left', bbox_to_anchor=(0.27, 1.12),fontsize=12)
AX[j].grid(zorder=0,ls=':',color='gray')
AX[j].set_yscale('log')
AX[j].set_xscale('log')
fig.supylabel(r'$N_{Halo(M_Stars)}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M))^{-1}$',fontsize=18)
fig.supxlabel('Log (Stellar Mass [M$_{\odot}$])',fontsize=18)
print()
plt.show()
print('\n\n\n')
fig= plt.figure(figsize=(17,5))
gs = gridspec.GridSpec( 1,5, # Number of axis y,x
height_ratios = [1], # relatives ratios of heigh
width_ratios = [1,0.1,1,0.1,1], # relatives ratio of with
left = 0.08, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.18, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2],sharey=ax1,sharex=ax1)
ax3 = fig.add_subplot(gs[0,4],sharey=ax1,sharex=ax1)
AX=[ax1,ax2,ax3]
# AX[0].ticklabel_format(axis='y', scilimits=[-1, 1])
AX[0].set_ylim(1e-5,1)
# Ploting the masses
for j,i in enumerate(LLL):
VV = (i)**3
XX = np.log10(globals()['mass_st_'+str(i)])+10##np.append(log_mass_st_q,log_mass_st_sf)
XX = XX[(XX>l)&(XX<u)]
y,x = np.histogram(XX,bins=50)
w = x[1]-x[0]
AX[j].plot(10**mid(x),y/VV/w,label=f'TNG {i}-1 (Stars)',marker='s',zorder=10)
AX[j].plot(10**F1,F2,label='SMP (Halo)',marker='.')
AX[j].plot(10**F3,F4,label='TNG 300 (Halo)',marker='.')
AX[j].plot(10**F5,F6,label='EAGLE (Halo)',marker='.')
AX[j].legend(loc='upper right',fontsize=12)#, bbox_to_anchor=(0.52, 1.12))
AX[j].grid(zorder=0,ls=':',color='gray')
AX[j].set_yscale('log')
AX[j].set_xscale('log')
fig.supylabel(r'$N_{Halo/Stars}\cdot\ dv^{-1}\cdot\ d(Log_{10}(M))^{-1}$',fontsize=18)
fig.supxlabel('Log (Mass [M$_{\odot}$])',fontsize=18)
print()
print('\n\n\n')
np.random.seed(1209)
Save_Galaxies={}
sub_data={}
for i in Data.keys():
# starformin and quenched subhalos [https://arxiv.org/pdf/1812.07584.pdf]
mass_st = globals()['mass_st_'+str(i)]
MH = Data[i]['mh']
sfr_tot = Data[i]['sfr_tot']
log_ssfr_tot = np.log10(sfr_tot / mass_st / 10**10)
mask_sf = (log_ssfr_tot > -11)
mask_q = (log_ssfr_tot < -11)
sub_data[i]=dict(log_mass_st_sf = np.log10(mass_st[mask_sf]) + 10,
log_mass_st_q = np.log10(mass_st[mask_q]) + 10,
log_sfr_tot_sf = np.log10(sfr_tot[mask_sf]),
log_sfr_tot_q = np.log10(sfr_tot[mask_q]),
MH_sf = MH[mask_sf],
MH_q = MH[mask_q])
# This is for activity 2.1
ids = Data[i]['sh_ids']
mask_mass = (-1<np.log10(mass_st))&(np.log10(mass_st)<1)
SF = np.random.choice(ids[(mask_sf & mask_mass)],10)
Q = np.random.choice(ids[(mask_q & mask_mass)],10)
Save_Galaxies[str(i)]={'SF':SF,'Q':Q}
if PRINT:
print(f'TNG {i}\nStar Forming: \t{SF}\nQuenched: \t{Q}')
print()
TNG 50 Star Forming: [632978 720707 386271 419620 627888 708395 167435 733149 300910 562743] Quenched: [504559 46 282781 470348 63937 96776 117259 289398 63875 143904] TNG 100 Star Forming: [577306 533722 419303 630829 595123 606391 580474 570763 332070 337846] Quenched: [134867 76208 279723 496713 182005 52859 114435 458028 197149 366793] TNG 300 Star Forming: [2116125 2168130 2511834 2256571 895130 1167096 1940368 1741389 2302679 1249631] Quenched: [ 902307 92281 834988 1232729 786460 904851 264577 777808 299342 1425898]
fig= plt.figure(figsize=(17,5))
gs = gridspec.GridSpec( 1,5, # Number of axis y,x
height_ratios = [1], # relatives ratios of heigh
width_ratios = [1,0.1,1,0.1,1], # relatives ratio of with
left = 0.08, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.18, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2],sharey=ax1,sharex=ax1)
ax3 = fig.add_subplot(gs[0,4],sharey=ax1,sharex=ax1)
AX=[ax1,ax2,ax3]
# AX[0].ticklabel_format(axis='y', scilimits=[-1, 1])
#AX[0].set_ylim(1e-5,1)
Var= ["log_mass_st_sf" ,"log_mass_st_q" , "log_sfr_tot_sf" ,"log_sfr_tot_q","MH_sf","MH_q"]
# Ploting the masses
for j,i in enumerate(LLL):
DDD = sub_data[str(i)]
AX[j].scatter(DDD[Var[0]], DDD[Var[2]], s=1, color='c', label='Star-Forming')
AX[j].scatter(DDD[Var[1]], DDD[Var[3]], s=1, color='m', label='Quenched')
AX[j].legend(loc='upper left',title=f'TNG{i}-1',fontsize=12, bbox_to_anchor=(.0, 1.12),framealpha=1)
AX[j].grid(zorder=0,ls=':',color='gray')
# AX[j].set_yscale('log')
# AX[j].set_xscale('log')
fig.supylabel('Log (SFR [M$_{\odot}$/yr])',fontsize=18)
fig.supxlabel('Log (Stellar Mass [M$_{\odot}$])',fontsize=18)
print()
plt.show()
print('\n\n\n')
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,3, # Number of axis y,x
height_ratios = [0.5,0.04,0.5], # relatives ratios of heigh
width_ratios = [0.33,0.33,0.33], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.1, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.2, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[2,0])
ax3 = fig.add_subplot(gs[0,1],sharex=ax1,sharey=ax1)
ax4 = fig.add_subplot(gs[2,1],sharex=ax2,sharey=ax2)
ax5 = fig.add_subplot(gs[0,2],sharex=ax1,sharey=ax1)
ax6 = fig.add_subplot(gs[2,2],sharex=ax2,sharey=ax2)
# (here we can define if we want to share axes and with which one)
AX=[ax1,ax2,ax3,ax4,ax5,ax6]
NNN= 9
#--
# Ploting
Var2=["sh_ids","cm","cv","rh_type","sfr_2rh","sfr_tot","mh"]
if PRINT:
print('#--------------------------------------------------------------')
print('#SET\tMedian Mass\t Median HMR\t MAD mass\t MAD HMR')
print('# \t Log(mass) \t Log(HMR) \tLog(mass)\t Log(HMR)')
for j,i in enumerate(LLL):
j=j*2
mass_st = np.log10(globals()['mass_st_'+str(i)])+10
HMR = np.log10(Data[str(i)][Var2[3]])
u,l= max(mass_st),min(mass_st)
AX[j].scatter(mass_st, HMR, s=1, color='b')#, label='HMR')
# AX[j].legend(loc='upper left',fontsize=12, bbox_to_anchor=(.0, .98),framealpha=1)
Lsp = np.linspace(l,u+0.1,NNN+1)
for O in range(NNN):
mask = (Lsp[O]<=mass_st)&(mass_st<Lsp[O+1])
x,y = np.median(mass_st[mask]),np.median(HMR[mask])
XE,YE = MAD(mass_st[mask]),MAD(HMR[mask])
if PRINT:
print(f'TNG{i:^4}: {x:.4f} \t {y:.4f} \t {XE:.4f} \t {YE:.4f}')
AX[j].errorbar(x,y, yerr=YE, xerr=XE,color='r')
DDD = sub_data[str(i)]
AX[j+1].scatter(DDD[Var[0]], np.log10(DDD[Var[4]]), s=1, color='c', label='Star Forming',alpha=0.5)
AX[j+1].scatter(DDD[Var[1]], np.log10(DDD[Var[5]]), s=1, color='m', label='Quenched',alpha=0.5)
AX[j].set_title(f'TNG{i}-1')
AX[j+1].legend(loc='lower right',fontsize=12, bbox_to_anchor=(1, 0),framealpha=1)
AX[j].grid(zorder=0,ls=':',color='gray')
AX[j+1].grid(zorder=0,ls=':',color='gray')
# AX[j].set_yscale('log')
# AX[j].set_xscale('log')
if PRINT:
print()
AX[1].set_ylim(-3.2,-1)
# fig.supylabel('Log (SFR [M$_{\odot}$/yr])',fontsize=18)
fig.supxlabel('Log (Stellar Mass [M$_{\odot}$])',fontsize=18)
# Press and Schechter
AX[0].set_ylabel("half-mass stellar radius Log(cKpc/h)",labelpad=12)
AX[1].set_ylabel("Metallicity Log(Z)",labelpad=12)
# Erasing extra y labels
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
plt.setp(ax5.get_xticklabels(), visible=False)
print()
plt.show()
print('\n\n\n')
#-------------------------------------------------------------- #SET Median Mass Median HMR MAD mass MAD HMR # Log(mass) Log(HMR) Log(mass) Log(HMR) TNG 50 : 7.9083 -0.1485 0.1402 0.1143 TNG 50 : 8.4990 -0.0979 0.1400 0.1451 TNG 50 : 9.0477 0.0990 0.1330 0.1645 TNG 50 : 9.6165 0.3084 0.1446 0.1763 TNG 50 : 10.2156 0.4381 0.1384 0.2115 TNG 50 : 10.7552 0.5779 0.1398 0.1613 TNG 50 : 11.2716 0.7661 0.1269 0.1085 TNG 50 : 11.8787 0.9771 0.1021 0.0907 TNG 50 : 12.7353 1.4666 0.0000 0.0000 TNG100 : 7.9000 0.3181 0.1392 0.0949 TNG100 : 8.4642 0.2740 0.1369 0.0852 TNG100 : 9.0157 0.2729 0.1358 0.0914 TNG100 : 9.5586 0.3215 0.1336 0.0980 TNG100 : 10.1140 0.3576 0.1354 0.1250 TNG100 : 10.6488 0.4736 0.1202 0.1585 TNG100 : 11.1681 0.8514 0.1179 0.1022 TNG100 : 11.7347 1.1998 0.1123 0.0862 TNG100 : 12.2431 1.5139 0.1043 0.1265 TNG300 : 7.8887 0.6524 0.1417 0.1237 TNG300 : 8.5125 0.5498 0.1491 0.1060 TNG300 : 9.1031 0.4868 0.1455 0.0839 TNG300 : 9.7022 0.4821 0.1470 0.0857 TNG300 : 10.2998 0.4918 0.1494 0.0781 TNG300 : 10.7860 0.6862 0.1098 0.1029 TNG300 : 11.3803 1.0967 0.1124 0.1074 TNG300 : 11.9668 1.5120 0.0949 0.1052 TNG300 : 12.5110 1.8852 0.0867 0.1512
if PRINT:
print(Save_Galaxies['50'])
{'SF': array([632978, 720707, 386271, 419620, 627888, 708395, 167435, 733149,
300910, 562743]), 'Q': array([504559, 46, 282781, 470348, 63937, 96776, 117259, 289398,
63875, 143904])}
# Some relevant functions
def get_basic_info_of_subhalo(Sim='50', ID=None):
DDD = Data[Sim]
if (Sim==None or ID==None):
print('please enter the index and ID of subhalo')
mask= np.where(DDD['sh_ids']==ID)[0][0]
subhalo_id = ID
r_half = DDD['rh_type'][mask] # half stellar mass radius
com = DDD['cm'][mask]
cov = DDD['cv'][mask]
return subhalo_id, r_half, com, cov
def cal_specific_ang_mom(mass, pos, vel, r_out):
r = LA.norm(pos, axis=1)
mask = np.where(r <= r_out)[0]
m = mass[mask]
p = pos[mask]
v = vel[mask]
L = np.cross(p, v)
L_total = np.sum(m * L.T, axis=1)
return L_total
def vector_rotation_matrix(vec1, vec2):
# unit vectors
v1 = vec1 / LA.norm(vec1)
v2 = vec2 / LA.norm(vec2)
# dimension and identity
dim = v1.size
I = np.identity(dim)
# cos of angle between v1 and v2
c = np.dot(v1, v2)
# cross product matrix of a vector to rotate around
K = np.outer(v2, v1) - np.outer(v1, v2)
# Rodrigues' formula
return I + K + (K @ K) / (1 + c)
i=50
basePath = '/home/tnguser/sims.TNG/TNG'+str(i)+'-1/output/'
Quen_pos={}
Quen_vel={}
Quen_posG={}
Quen_velG={}
DM_or_gas_Q={}
MASS_DM = ((3.64755609660833e-05)/h)*(1e10) # https://www.tng-project.org/data/downloads/TNG50-1-Dark/
if PRINT:
print(MASS_DM)
for i in Save_Galaxies['50']['Q']:
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,7, # Number of axis y,x
height_ratios = [1,0.23,1], # relatives ratios of heigh
width_ratios = [1,0.1,1,0.1,1,0.1,0.08], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.15, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2])
ax3 = fig.add_subplot(gs[0,4])
cax1 = fig.add_subplot(gs[0,6])
ax4 = fig.add_subplot(gs[2,0])
ax5 = fig.add_subplot(gs[2,2])
ax6 = fig.add_subplot(gs[2,4])
cax2 = fig.add_subplot(gs[2,6])
axs=[ax1,ax2,ax3,ax4,ax5,ax6]
fig.suptitle(f'Galaxy {i} / Quenched',y=1.05)
#################################################################################
subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(ID=i)
fields = ['Coordinates', 'Velocities', 'Masses']
st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'stars',fields=fields)
mass = st['Masses']
pos = st['Coordinates']
vel = st['Velocities']
pos = pos - com # shifting co-ordinate system to center of subhalo
vel = vel - cov
# plotting three projection of subhalo after rotation
#fist rotate subhalo
L_vec = cal_specific_ang_mom(mass, pos, vel, 10*r_half) # angular momentum vector
rot_mat = vector_rotation_matrix(L_vec, [0,0,1]) # [0,0,1] for face-on
new_pos = np.dot(rot_mat, pos.T)
new_vel = np.dot(rot_mat, vel.T)
new_pos, new_vel = new_pos.T, new_vel.T
Quen_pos[str(i)]=new_pos
Quen_vel[str(i)]=new_vel
extent = 5 # in unit of half mass radius
rng = [[-extent*r_half,extent*r_half],[-extent*r_half,extent*r_half]]
nbins = int((rng[0][1] - rng[0][0])/0.12)
H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,1], bins=nbins, range=rng)
H1[H1<=0]=np.nan
axs[0].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[0].set_xlabel('x [kpc]')
axs[0].set_ylabel('y [kpc]')
H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,2], bins=nbins, range=rng)
H1[H1<=0]=np.nan
axs[1].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[1].set_xlabel('x [kpc]')
axs[1].set_ylabel('z [kpc]')
H1,X1,Y1= np.histogram2d(new_pos[:,1], new_pos[:,2], bins=nbins, range=rng)
H1[H1<=0]=np.nan
IM=axs[2].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[2].set_xlabel('y [kpc]')
axs[2].set_ylabel('z [kpc]')
plt.colorbar(IM,cax=cax1,ticks=[np.log10(3)*1.05,np.log10(np.nanmax(H1))*0.95])
cax1.set_ylabel('Density', rotation=270)
cax1.set_title('Stars')
cax1.set_yticklabels(['low','high'])
####################################################################################
####################################################################################
############## GAS ###########################################################
####################################################################################
####################################################################################
try:
st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'gas',fields=fields)
if st['count']<5000:
raise Exception('To few particles showing instead DM')
mass = st['Masses']
pos = st['Coordinates']
vel = st['Velocities']
cax2.set_title('Gas')
print('Gas: Count: ',st['count'])
CMAP = 'inferno'
DM_or_gas_Q[str(i)]='Gas'
except Exception as e:
print(e)
print('Gas: Count: ',st['count'])
fields = ['Coordinates', 'Velocities']
st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'dm',fields=fields)
pos = st['Coordinates']
vel = st['Velocities']
mass = np.repeat(MASS_DM,len(pos))
cax2.set_title('DM')
CMAP = 'viridis'
DM_or_gas_Q[str(i)]='Dark Matter'
pos = pos - com # shifting co-ordinate system to center of subhalo
vel = vel - cov
# plotting three projection of subhalo after rotation
#fist rotate subhalo
L_vec = cal_specific_ang_mom(mass, pos, vel, 10*r_half) # angular momentum vector
rot_mat = vector_rotation_matrix(L_vec, [0,0,1]) # [0,0,1] for face-on
new_pos = np.dot(rot_mat, pos.T)
new_vel = np.dot(rot_mat, vel.T)
new_pos, new_vel = new_pos.T, new_vel.T
Quen_posG[str(i)]=new_pos
Quen_velG[str(i)]=new_vel
H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,1], bins=nbins, range=rng)
H1[H1<=0]=np.nan
axs[3].imshow(np.log10(H1.T) ,vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[3].set_xlabel('x [kpc]')
axs[3].set_ylabel('y [kpc]')
H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,2], bins=nbins, range=rng)
H1[H1<=0]=np.nan
axs[4].imshow(np.log10(H1.T) ,vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[4].set_xlabel('x [kpc]')
axs[4].set_ylabel('z [kpc]')
H1,X1,Y1= np.histogram2d(new_pos[:,1], new_pos[:,2], bins=nbins, range=rng)
H1[H1<=0]=np.nan
IM=axs[5].imshow(np.log10(H1.T),vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[5].set_xlabel('y [kpc]')
axs[5].set_ylabel('z [kpc]')
plt.colorbar(IM,cax=cax2,ticks=[np.log10(np.nanmin(H1))*1.07,np.log10(np.nanmax(H1))*0.95])
cax2.set_ylabel('Density', rotation=270)
cax2.set_yticklabels(['low','high'])
plt.show()
print('\n')
print('\n\n\n')
538464.1418081385 Gas: Count: 494611
To few particles showing instead DM Gas: Count: 0
To few particles showing instead DM Gas: Count: 4102
To few particles showing instead DM Gas: Count: 0
To few particles showing instead DM Gas: Count: 0
To few particles showing instead DM Gas: Count: 7
To few particles showing instead DM Gas: Count: 486
To few particles showing instead DM Gas: Count: 0
To few particles showing instead DM Gas: Count: 0
To few particles showing instead DM Gas: Count: 0
SFR_pos={}
SFR_vel={}
SFR_posG={}
SFR_velG={}
DM_or_gas_SF={}
for i in Save_Galaxies['50']['SF']:
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,7, # Number of axis y,x
height_ratios = [1,0.23,1], # relatives ratios of heigh
width_ratios = [1,0.1,1,0.1,1,0.1,0.08], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.15, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2])
ax3 = fig.add_subplot(gs[0,4])
cax1 = fig.add_subplot(gs[0,6])
ax4 = fig.add_subplot(gs[2,0])
ax5 = fig.add_subplot(gs[2,2])
ax6 = fig.add_subplot(gs[2,4])
cax2 = fig.add_subplot(gs[2,6])
axs=[ax1,ax2,ax3,ax4,ax5,ax6]
fig.suptitle(f'Galaxy {i} / Star Forming',y=1.05)
#################################################################################
subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(ID=i)
fields = ['Coordinates', 'Velocities', 'Masses']
st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'stars',fields=fields)
mass = st['Masses']
pos = st['Coordinates']
vel = st['Velocities']
pos = pos - com # shifting co-ordinate system to center of subhalo
vel = vel - cov
# plotting three projection of subhalo after rotation
#fist rotate subhalo
L_vec = cal_specific_ang_mom(mass, pos, vel, 10*r_half) # angular momentum vector
rot_mat = vector_rotation_matrix(L_vec, [0,0,1]) # [0,0,1] for face-on
new_pos = np.dot(rot_mat, pos.T)
new_vel = np.dot(rot_mat, vel.T)
new_pos, new_vel = new_pos.T, new_vel.T
SFR_pos[str(i)]=new_pos
SFR_vel[str(i)]=new_vel
extent = 5 # in unit of half mass radius
rng = [[-extent*r_half,extent*r_half],[-extent*r_half,extent*r_half]]
nbins = int((rng[0][1] - rng[0][0])/0.12)
H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,1], bins=nbins, range=rng)
H1[H1<=0]=np.nan
axs[0].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[0].set_xlabel('x [kpc]')
axs[0].set_ylabel('y [kpc]')
H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,2], bins=nbins, range=rng)
H1[H1<=0]=np.nan
axs[1].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[1].set_xlabel('x [kpc]')
axs[1].set_ylabel('z [kpc]')
H1,X1,Y1= np.histogram2d(new_pos[:,1], new_pos[:,2], bins=nbins, range=rng)
H1[H1<=0]=np.nan
IM=axs[2].imshow(np.log10(H1.T), vmin=np.log10(3) ,cmap='gist_ncar',interpolation = 'spline36',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[2].set_xlabel('y [kpc]')
axs[2].set_ylabel('z [kpc]')
plt.colorbar(IM,cax=cax1,ticks=[np.log10(3)*1.05,np.log10(np.nanmax(H1))*0.95])
cax1.set_ylabel('Density', rotation=270)
cax1.set_title('Stars')
cax1.set_yticklabels(['low','high'])
####################################################################################
####################################################################################
############## GAS or DM ######################################################
####################################################################################
####################################################################################
try:
st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'gas',fields=fields)
if st['count']<100:
raise Exception('To few particles showing instead DM')
mass = st['Masses']
pos = st['Coordinates']
vel = st['Velocities']
cax2.set_title('Gas')
print('Gas: Count: ',st['count'])
CMAP = 'inferno'
DM_or_gas_SF[str(i)]='Gas'
except Exception as e:
print(e)
print('Gas: Count: ',st['count'])
fields = ['Coordinates', 'Velocities']
st = il.snapshot.loadSubhalo(basePath,snapNum,subhalo_id,'dm',fields=fields)
pos = st['Coordinates']
vel = st['Velocities']
mass = np.repeat(MASS_DM,len(pos))
cax2.set_title('DM')
CMAP = 'viridis'
DM_or_gas_SF[str(i)]='Dark Matter'
pos = pos - com # shifting co-ordinate system to center of subhalo
vel = vel - cov
# plotting three projection of subhalo after rotation
#fist rotate subhalo
L_vec = cal_specific_ang_mom(mass, pos, vel, 10*r_half) # angular momentum vector
rot_mat = vector_rotation_matrix(L_vec, [0,0,1]) # [0,0,1] for face-on
new_pos = np.dot(rot_mat, pos.T)
new_vel = np.dot(rot_mat, vel.T)
new_pos, new_vel = new_pos.T, new_vel.T
SFR_posG[str(i)]=new_pos
SFR_velG[str(i)]=new_vel
H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,1], bins=nbins, range=rng)
H1[H1<=0]=np.nan
axs[3].imshow(np.log10(H1.T) ,vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[3].set_xlabel('x [kpc]')
axs[3].set_ylabel('y [kpc]')
H1,X1,Y1= np.histogram2d(new_pos[:,0], new_pos[:,2], bins=nbins, range=rng)
H1[H1<=0]=np.nan
axs[4].imshow(np.log10(H1.T) ,vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[4].set_xlabel('x [kpc]')
axs[4].set_ylabel('z [kpc]')
H1,X1,Y1= np.histogram2d(new_pos[:,1], new_pos[:,2], bins=nbins, range=rng)
H1[H1<=0]=np.nan
IM=axs[5].imshow(np.log10(H1.T),vmin=np.log10(np.nanmin(H1))*1.05,cmap=CMAP,interpolation = 'spline16',interpolation_stage='rgba',origin='lower',extent=[-extent*r_half,extent*r_half,-extent*r_half,extent*r_half])
axs[5].set_xlabel('y [kpc]')
axs[5].set_ylabel('z [kpc]')
plt.colorbar(IM,cax=cax2,ticks=[np.log10(np.nanmin(H1))*1.07,np.log10(np.nanmax(H1))*0.95])
cax2.set_ylabel('Density', rotation=270)
cax2.set_yticklabels(['low','high'])
plt.show()
print('\n')
print('\n\n\n')
Gas: Count: 351436
Gas: Count: 71575
Gas: Count: 349190
Gas: Count: 110659
Gas: Count: 358181
Gas: Count: 67241
Gas: Count: 15699
Gas: Count: 65494
Gas: Count: 5464
Gas: Count: 126201
## To save the data of the galaxies
# np.save('SFR_pos.npy', SFR_pos)
# np.save('SFR_vel.npy', SFR_vel)
# np.save('Quen_pos.npy', Quen_pos)
# np.save('Quen_vel.npy', Quen_vel)
N=100
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,9, # Number of axis y,x
height_ratios = [1,0.2,1], # relatives ratios of heigh
width_ratios = [1,0.05,1,0.05,1,0.05,1,0.05,1], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.15, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2],sharey=ax1)
ax3 = fig.add_subplot(gs[0,4],sharey=ax1)
ax4 = fig.add_subplot(gs[0,6],sharey=ax1)
ax5 = fig.add_subplot(gs[0,8],sharey=ax1)
ax6 = fig.add_subplot(gs[2,0],sharey=ax1)
ax7 = fig.add_subplot(gs[2,2],sharey=ax1)
ax8 = fig.add_subplot(gs[2,4],sharey=ax1)
ax9 = fig.add_subplot(gs[2,6],sharey=ax1)
ax10 = fig.add_subplot(gs[2,8],sharey=ax1)
AX=[ax1,ax2,ax3,ax4,ax5,ax6,ax7,ax8,ax9,ax10]
AX[0].set_yscale('log')
COLOR={'Stars':'blue','Dark Matter':'springgreen','Gas':'red'}
for j,i in enumerate(Quen_pos.keys()):
subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(ID=int(i))
#####################################################
# Stars
#####################################################
Dato = Quen_pos[i]
R_part= (Dato[:,0]**2+Dato[:,1]**2+Dato[:,2]**2)**0.5
R_part= R_part[(R_part<=r_half)]
H,E = np.histogram(R_part,bins=N)
W= ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3))
#V= (4/3)*np.pi*r_half**3
AX[j].stairs(H/W,E/r_half,label='Stars',color=COLOR['Stars']) # Divided by the volume enclosed , normilized to the half radius
AX[j].set_title(f'Galaxy {i}/ Quenched',fontsize=10)
#####################################################
# Gas or DM
#####################################################
Dato = Quen_posG[i]
TYPE = DM_or_gas_Q[str(i)]
R_part= (Dato[:,0]**2+Dato[:,1]**2+Dato[:,2]**2)**0.5
R_part= R_part[(R_part<=r_half)]
H,E = np.histogram(R_part,bins=N)
W= ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3))
#V= (4/3)*np.pi*r_half**3
AX[j].stairs(H/W,E/r_half,label=TYPE,color=COLOR[TYPE]) # Divided by the volume enclosed , normilized to the half radius
LG=AX[j].legend(loc='upper right',fontsize=9)
LG.set_title(f'$R_{{half}}: {r_half:.3f}$',prop={'size':10})
# Erasing extra y labels
plt.setp(ax2.get_yticklabels(), visible=False)
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)
plt.setp(ax5.get_yticklabels(), visible=False)
plt.setp(ax7.get_yticklabels(), visible=False)
plt.setp(ax8.get_yticklabels(), visible=False)
plt.setp(ax9.get_yticklabels(), visible=False)
plt.setp(ax10.get_yticklabels(), visible=False)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
plt.setp(ax4.get_xticklabels(), visible=False)
plt.setp(ax5.get_xticklabels(), visible=False)
fig.supylabel('Log ([$N_{particules}/kpc^3$])',fontsize=18)
fig.supxlabel('Log ($R/R_{Half}$)',fontsize=18)
plt.show()
print('\n')
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 3,9, # Number of axis y,x
height_ratios = [1,0.2,1], # relatives ratios of heigh
width_ratios = [1,0.05,1,0.05,1,0.05,1,0.05,1], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.15, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2],sharey=ax1)
ax3 = fig.add_subplot(gs[0,4],sharey=ax1)
ax4 = fig.add_subplot(gs[0,6],sharey=ax1)
ax5 = fig.add_subplot(gs[0,8],sharey=ax1)
ax6 = fig.add_subplot(gs[2,0],sharey=ax1)
ax7 = fig.add_subplot(gs[2,2],sharey=ax1)
ax8 = fig.add_subplot(gs[2,4],sharey=ax1)
ax9 = fig.add_subplot(gs[2,6],sharey=ax1)
ax10 = fig.add_subplot(gs[2,8],sharey=ax1)
AX=[ax1,ax2,ax3,ax4,ax5,ax6,ax7,ax8,ax9,ax10]
AX[0].set_yscale('log')
for j,i in enumerate(SFR_pos.keys()):
subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(ID=int(i))
#####################################################
# Stars
#####################################################
Dato = SFR_pos[i]
R_part= (Dato[:,0]**2+Dato[:,1]**2+Dato[:,2]**2)**0.5
R_part= R_part[(R_part<=r_half)]
H,E = np.histogram(R_part,bins=N)
W= ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3))
#V= (4/3)*np.pi*r_half**3
AX[j].stairs(H/W,E/r_half,label='Stars',color=COLOR['Stars'])
AX[j].set_title(f'Galaxy {i}/ Star Forming',fontsize=10)
#####################################################
# Gas or DM
#####################################################
Dato = SFR_posG[i]
TYPE = DM_or_gas_SF[str(i)]
R_part= (Dato[:,0]**2+Dato[:,1]**2+Dato[:,2]**2)**0.5
R_part= R_part[(R_part<=r_half)]
H,E = np.histogram(R_part,bins=N)
W= ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3))
#V= (4/3)*np.pi*r_half**3
AX[j].stairs(H/W,E/r_half,label=TYPE,color=COLOR[TYPE]) # Divided by the volume enclosed , normilized to the half radius
LG=AX[j].legend(loc='upper right',fontsize=9)
LG.set_title(f'$R_{{half}}: {r_half:.3f}$',prop={'size':10})
# Erasing extra y labels
plt.setp(ax2.get_yticklabels(), visible=False)
plt.setp(ax3.get_yticklabels(), visible=False)
plt.setp(ax4.get_yticklabels(), visible=False)
plt.setp(ax5.get_yticklabels(), visible=False)
plt.setp(ax7.get_yticklabels(), visible=False)
plt.setp(ax8.get_yticklabels(), visible=False)
plt.setp(ax9.get_yticklabels(), visible=False)
plt.setp(ax10.get_yticklabels(), visible=False)
plt.setp(ax1.get_xticklabels(), visible=False)
plt.setp(ax2.get_xticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)
plt.setp(ax4.get_xticklabels(), visible=False)
plt.setp(ax5.get_xticklabels(), visible=False)
fig.supylabel('Log ([$N_{particules}/kpc^3$])',fontsize=18)
fig.supxlabel('Log ($R/R_{Half}$)',fontsize=18)
plt.show()
print('\n')
# import astropy.units as u
# p_crit=((1.9e-29)*(h**2))*u.g*(u.cm**-3)
# p_crit=p_crit.to(u.Msun*(u.kpc**-3))
# p_crit
# p_crit=p_crit.value
i = 100
basePath = '/home/tnguser/sims.TNG/TNG'+str(i)+'-1/output/'
FIELDS=['GroupFirstSub','GroupMass','GroupMassType','Group_M_Crit200','Group_R_Crit200']
HHH=il.groupcat.loadHalos(basePath,99,fields=FIELDS)
GFSub = HHH['GroupFirstSub']
GMass = HHH['GroupMass']
Mtype = HHH['GroupMassType']
M200 = HHH['Group_M_Crit200']
R200 = HHH['Group_R_Crit200']
# https://www.tng-project.org/data/downloads/TNG100-1/
h = 0.6774
MASS_DM = (0.000505574296436975e10)/h # M sun
Halos = {}
ddd = GMass.copy()
if PRINT:
print('# Most Massive v/s Group First Sub')
for i in range(3):
a=ddd.argmax()
if PRINT:
print(f' {a}\t{GFSub[i]}')
ddd[a]=0
#pos = il.groupcat.loadSingle(basePath,99,subhaloID=GroupFirstSub[i])
subhalo_id, r_half, com, cov = get_basic_info_of_subhalo(Sim='100',ID=a)
pos = il.snapshot.loadSubhalo(basePath,99,a,'dm',fields=['Coordinates'])
pos = pos-com
R_part = ((pos[:,0]**2+pos[:,1]**2+pos[:,2]**2)**0.5)/h
Halos[a] = R_part
del ddd
if PRINT:
print(Halos.keys())
# Most Massive v/s Group First Sub 0 0 1 17185 2 31342 dict_keys([0, 1, 2])
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 1,5, # Number of axis y,x
height_ratios = [1], # relatives ratios of heigh
width_ratios = [1,0.05,1,0.05,1], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.15, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2])#,sharey=ax1)
ax3 = fig.add_subplot(gs[0,4])#,sharey=ax1)
AX=[ax1,ax2,ax3]
N=100
XY_data={}
for j,i in enumerate(Halos.keys()):
R_part = Halos[i] # Get the data stored
R_trim = R_part[((0.74<R_part)&(R_part<R200[i]))] # Cut the data in the range asked
R_trim = np.log10(R_trim) # in log space
H,E = np.histogram(R_trim,bins=N)
E = 10**E
W = ((4/3)*np.pi*(E[1:]**3)) - ((4/3)*np.pi*(E[:-1]**3)) # Volume in each step
Y = (H/W)*MASS_DM
XY_data[i]=[mid(E),Y]
AX[j].stairs(Y,E/R200[i]
,label=f'Halo {i}')
#AX[j].axhline(200*p_crit,color='red')
AX[j].set_yscale('log')
AX[j].set_xscale('log')
AX[j].legend(title=f'$R_{{200}}$= {R200[i]:.2f}',loc='upper right')
fig.supylabel('Log ([M$_{\odot}/kpc^3$])',fontsize=18)
fig.supxlabel('Log (R/$R_{200}$)',fontsize=18)
plt.show()
print('\n\n\n')
# We will use the scipy package in order to get the best fit
from scipy.optimize import curve_fit
def NFW(x,p,r): # Navarro-Frenk-White (NFW) profile.
X=(x/r)
return np.log10(p/(X*((1+X)**2))) # We do the fit in Log space
fig= plt.figure(figsize=(18,9))
gs = gridspec.GridSpec( 1,5, # Number of axis y,x
height_ratios = [1], # relatives ratios of heigh
width_ratios = [1,0.05,1,0.05,1], # relatives ratio of with
left = 0.1, # Space to the edge of left from the nearest axis
right = 0.95, # Space to the edge of right from the nearest axis
bottom= 0.15, # Space to the edge of bottom from the nearest axis
top = 0.97, # Space to the edge of top from the nearest axis
wspace= 0.1, # Space horizontal between each of the axis
hspace= 0.0) # Space vertical between each of the axis
ax1 = fig.add_subplot(gs[0,0])
ax2 = fig.add_subplot(gs[0,2])#,sharey=ax1)
ax3 = fig.add_subplot(gs[0,4])#,sharey=ax1)
AX=[ax1,ax2,ax3]
AAA=[[1e6,300],[1e6,150],[1e7,40]] #Init guess (something reasonable)
for j,i in enumerate(XY_data.keys()):
xdata=(XY_data[i][0])
ydata=np.log10(XY_data[i][1]) # We do the fit in Log space
A=AAA[j]
popt, pcov = curve_fit(NFW, xdata, ydata,maxfev=10000,p0=A,bounds=(0.74, [1e16, 1e3]))
if PRINT:
print(f'For Halo {i} we have a $rho_0$= {popt[0]:.3f},\t and a $R_s$={popt[1]:.3f}')
AX[j].scatter(xdata/R200[i] ,ydata -np.log10(popt[0]) ,label='Data', marker='o',facecolor='none',edgecolor='r')
AX[j].plot( xdata/R200[i] ,NFW(xdata,popt[0],popt[1])-np.log10(popt[0]) ,label='Best Fit')
AX[j].plot( xdata/R200[i] ,NFW(xdata,A[0],A[1]) -np.log10(popt[0]) ,label='My eye Fit')
# plt.yscale('log')
AX[j].set_xscale('log')
AX[j].set_title(f'Halo {i}')
LG=AX[j].legend(loc='upper right',fontsize=9)
LG.set_title('From Best Fit\n'+rf'Log ($\rho_0$)= {np.log10(popt[0]):.3f}'+f'\n$R_s$ = {popt[1]:.2f}',prop={'size':10})
fig.supylabel(r'Log ($\rho/\rho_0$)',fontsize=18)
fig.supxlabel(r'Log ($R/R_0$)',fontsize=18)
print()
For Halo 0 we have a $rho_0$= 552491.892, and a $R_s$=365.993 For Halo 1 we have a $rho_0$= 4189014.063, and a $R_s$=69.198 For Halo 2 we have a $rho_0$= 48541611.218, and a $R_s$=14.726
print('\n\n\nBy: Ian Baeza')
By: Ian Baeza